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Vertical density matrix algorithm (VDMA), a tensor product state formulation of the "higher- 
dimensional" density matrix renormalization group, is applied to the spin 1/2 antiferromagnetic 
XXZ model on the checkerboard lattice. The VDMA was, in the preceding study, applied to 
the transverse field Ising model on the square lattice and the three-dimensional classical Ising 
model. In the present paper, its implementation procedure is modified in order to apply the 
VDMA to the XXZ model. Numerical accuracy of the VDMA is investigated for the XXZ 
model on the square lattice, which shows that the method gives reliable results for the ground 
state energy. In the frustrated region, VDMA results are compared with a simple calculation 
based on a magnetically disordered state. It is found that the weakly frustrated region is in 
the Neel ordered phase, while in the strongly frustrated region the realized phase cannot be 
identified clearly by the obtained results. 

KEYWORDS: DMRG, tensor product state, quantum spin system, frustration 



1. Introduction 

The density matrix renormalization group (DMRG) 
method^' ^ has become one of the standard numerical 
techniques for studying one-dimensional (ID) quantum 
systems and 2D classical systems because of numerical 
accuracy and the ability to deal with large size systems.'^ 
The success of the DMRG has been stimulating us to 
extend the algorithm to the one that can handle higher- 
dimensional systems, principally 2D quantum systems 
and 3D classical sy stems. Nevertheless, we have not 
achieved the definite success. 

One of the key ideas for developing the "higher- 
dimensional" DMRG is the tensor product state 
(TPS).^"^^ The TPS is a higher-dimensional generaliza- 
tion of the matrix product state (MPS) for ID quantum 
spin systems. On the basis of the idea of the TPS, sev- 
eral formulations for 3D classical systems have been pro- 
posed. "'^^"^^ One of these TPS formulations is the vertical 
density matrix algorithm (VDMA).^^ In the preceding 
study, it is demonstrated that, besides the 3D classical 
Ising model, the 2D transverse field Ising (TFI) model, a 
most fundamental 2D quantum spin system, can be also 
dealt with by using VDMA.^^ However, its implementa- 
tion procedure for the 2D TFI model is completely equiv- 
alent to the one for the 3D classical Ising model because 
the 2D TFI model is exactly mapped to the 3D classical 
Ising model. Hence further development and modifica- 
tions for general 2D quantum systems are desired for the 
progress of the higher-dimensional algorithm. 

In this paper, the VDMA is applied to the S=l/2 an- 
tiferromagnetic XXZ model on the checkerboard lattice. 
This lattice consists of the corner-sharing plaquettes that 
have the nearest neighbor and the next-nearest neighbor 
bonds, as depicted in Fig. 1. Interestingly, in the isotropic 
case, this model is predicted to exhibit either the Neel 
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Fig. 1. The checkerboard lattice. Solid (dashed) lines denote the 
nearest (next-nearest) neighbor bonds. 



ordered phase^^"^* or a magnetically disordered phase 
called the valence bond crystal (VBC).'^^"'^^ Hence it is 
appropriate to test the capability of the VDMA for both 
cases. 

First, the VDMA calculation has been performed for 
the unfrustrated 2D XXZ model to compare the VDMA 
results with available quantum Monte Carlo (QMC) 
data. It is found that the VDMA gives quantitatively 
reliable results for the ground state energy; the relative 
deviations from the QMC results are the order of 1%. In 
the frustrated region, the VDMA results are compared 
with a simple calculation based on the VBC. It is found 
that in the weakly frustrated region the Neel ordered 
phase exists, while in the strongly frustrated region the 
realized phase cannot be identified by the obtained re- 
sult; in the latter region, the VDMA result for the energy 
is greater than the upper bound of the energy obtained 
by using the VBC state. 

This paper is organized as follows: In § 2, the definition 
of the model is given, and then the VDMA for the 2D 
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Fig. 2. (a): the checkerboard lattice with plaquettes divided into 
two groups A and B. (b): the plaquette labelled siij). 
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fig. 3 N. Maeshima 

Fig. 3. The lattice structure of the 3D classical model mapped 
from the 2D XXZ model. The Z-axis corresponds to the Trotter 
axis. 

XXZ model on the checkerboard lattice is explained. Sec- 
tion 3 shows the numerical results for the unfrustrated 
case and the frustrated case. The last section is devoted 
to the conclusion. 

2. The Vertical Density Matrix Algorithm for 
the 2D XXZ Model 

The outline of the VDMA for 2D quantum spin sys- 
tems is given in ref.l7; after the mapping of 2D quan- 
tum spin systems to 3D classical spin systems by the 
Suzuki- Trotter decomposition,"^^ the VDMA is used to 
the 3D classical systems. However, the implementation 
procedure for the 2D XXZ model is significantly different 
from the one for the 2D TFI model, which is described in 
ref.l7, because of the difference in the lattice structure 
of the mapped 3D classical systems. 

Let H be the Hamiltonian of the 2D quantum spin 
model under consideration; its explicit form is given as 

(W).(u)> 

The sum runs over the nearest bonds and the next- 
nearest bonds shown in Fig. 1, and is the xiy) coor- 
dinate of a site. The local Hamiltonian with the 
anisotropy parameter g'(>0) is defined as 

HaW) = 9{SijSfj + SfjSl.) + SljSf-., (2) 



where Sf^ is the a-component (a = x, y, or z) of the spin 
operator at the site ij. The exchange coupling J(ij)(ij) 
is set to unity for the nearest neighbor bonds and jd for 
the next nearest neighbor bonds. 

The first step of the algorithm is to map the 2D quan- 
tum spin model into a 3D classical spin model. We fol- 
low the standard "checkerboard" decomposition^" for 2D 
quantum systems. Divide all the corner-sharing plaque- 
ttes in the checkerboard lattice into two groups A and B, 
as shown in Fig. 2 (a). A plaquette that has 4 sites {{ij), 
{ij'), {i'j), and {i'j' j) is labeled as s{ij), where i' = i + 1 
and j' = j + 1 (see Fig. 2(b)). Then the Hamiltonians of 
the two groups A and B are written as 

Ha = ^ ^s{ij) ) and Hb = ^ ^s(ij) ■ (3) 

s{ij)£A s{ij)GB 

Here hg^ij-j is the local Hamiltonian of the 4 spins on the 
plaquette s{ij); 

+ 3d[h^3){i'r)+^(i'3m')\- (4) 
Then the following Suzuki- Trotter transformation 

Z= lim Tr(e-'^-^e-'^«)^ (5) 

is used to map the 2D XXZ model at temperature T into 
the classical model on the 3D lattice with 2L Trotter 
size (see Fig. 3). It is found that the mapped 3D lattice 
is composed of the corner-sharing cubes, each of which 
have 8 Ising spins at its corners. The Boltzmann weight 
for an unit cube is represented as 

where we use the notation e = J/{LkBT) and |(Ty ) — 
■ Wij) is the ^^-representation of the 

spin at the site ij. 

The next step is to take the L — »■ oo limit with fixed e, 
which is equivalent to taking the zero temperature limit 
of the original 2D model. Then the partition function Z 
can be expressed as a simple form. Consider e~^^^e~^^® 
as a matrix T, and let -Bmax and liJjmaL) be the maximum 
eigenvalue and the corresponding right (left) eigenstate 
of T. Using the power method, we can calculate the par- 
tition function as follows; 

Z = lim TrT^ 

Z/—J-00, fixed e 

= , ^ i^Lx(V'La.lCa.>- (7) 

-L— >c»,nxea e 

Hence it is found that Z is equivalent to the norm 
(V'maxlV'max) ^xcept for the trivial coefficient -E^ax- Here 
it should be stressed that we now need both l^/'^ax) ^^'^ 
(V'max I to evaluate the partition function because T is an 
asymmetric matrix. In contrast, for the 2D TFI model 
we can use a symmetric matrix, and therefore we needs 
only the right (or left) eigenstate.^'' 

In the framework of the VDMA, each eigenstate is 
given in a form of the TPS. Let us consider the TPS 
structure of IV'max) using the power method. Define a 
vector |'^£) as 

\rL)=T'^\i'o), (8) 
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where I^Ao) is an "initial" vector that is not orthogonal fig-4 N. Maeshima 
to IV'max)- Here we use given by a TPS; 



s(ij)G-B 



(9) 



where w(Tf^) is a local tensor. Then |'0£} is rewritten as 
the TPS; 



E 

n 

s(ii)e-B 



n ^ 



t2L 



(#) 



where 



2L\ 

ij I 



(10) 



(11) 



The two tensors and are defined as 

,2L 



and 



^2X 



W 



r3 



w 



2L-2 
,2L-3 



,1 



(12) 



j2L 
'ij 



w 



2L-1 

ij 



W 



W 



w{tI), (13) 



where $,f^ = {crl^ afj ■■■ cr'^^ ^) is an auxiliary vari- 
able." 

It should be noted that | ■)/'£) consists of two types of 
tensors X"^ and X^ whereas only one tensor is needed 
for the 2D TFI model. This difference stems from the 
structure of |'0£), which is shown in Fig. 4. We can see 
that It/i^) is constructed by X^ and X^ , which are rep- 
resented by rectangular parallelepipeds composed of the 
cubes vertically piling up (see Fig. 5). We also found that 
the bare spin variables {o'tj"} are on the upper faces of 
the rectangular parallelepipeds corresponding to the ten- 
sor X^, not to X^ . Hence the tensor X"^ has the bare 
spin variables afj" as indices and X^ does not so. 

By taking the limit L —^ oo, the right maximum- 
eigenvalue eigenvector li/^max) is obtained; 



^) = const X lim 



(14) 



Therefore IV'max) is represented as the TPS composed of 
the tensors X^ and X^ with the auxiliary variables . 
The left maximum-eigenvalue eigenvector (V'maxI is also 
represented as the same TPS form with the "transposed 
local tensor" , 



X' 



w 



w(t%)W 



w 



w 



- 2L-2 

ij 
-2L-1 



r3 



(15) 
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Fig. 4. The graphical representation of the vector l'/'^} {L=4). 
Filled circles and crosses show the summed spin variables in 
eq. (10). Open circles represent the fixed spin variables. 
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Fig. 5. The two tensors JsT^ and . 
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where ryff = (o-i-cr 



tor 



^-2 ) '''' I _iL ] ; (16) 

^). When we define the vec- 



(^Ll = (V'o|T^ (17) 
('/'max I is represented as the following TPS form; 

(V-Lxlk""]) = const X lim (^^^[(t^^]) (18) 



= const X 



hm y 

^.i-allij 



n ^l{v 

s{i])eA 



2L\ 
ij ) 



X 



n > 



2L 



(19) 



It should be also noted that, from the definitions of the 
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tensors, the following relations hold; 



X' 



(20) 



These relations are used to obtain renormalized tensors, 
as is mentioned later. 

The key point of the VDMA is to reduce the dimen- 
sions of the auxiliary variables ^ and rf using the vertical 
density matrix (VDM) p. By taking account of the lat- 
tice structure of the 3D classical model, we can define 
the VDM as 

Pkl{crkl£,'klWkl£.kl) = 



E' 

n 

s{if)eB 



X' 



X 



ik"l" 



X' 



Vkl 
CTkl 



X^iVk"i") 



X^'i^ki). (21) 



where ^ denotes the configuration sum for the all 
spin variables except o-^i^^ki^'^ki, and ^u- Yi 



s{ij)GA 



de- 



notes the product for the plaquettes of group A except 
s{k",l") = s{k — 1,1 — 1), and Y[ s(ij)GB product 
for the plaquettes of group B except s{k, I). The checked 
spin is defined as 



f fe"/' 

^k"l" 



Ck"l" 



Ckl'- 
^kl" 



'kl 
kl 



£,k"l 



(22) 



The locations of the spin variables are shown in Fig. 6 

The diagonalization of the VDM p gives us the pro- 
jection operators to reduce the dimensions of ^ and 
rj. Since the VDM p is asymmetric, we obtain right- 
eigenvectors [/'"(cr^lC) and the left-eigenvectors U\a£_\£,), 
with the corresponding eigenvalues uj^ in the decreasing 
order lui > lu2 > ■ • • . In the actual implementation of 
the VDMA for the 2D XXZ model, all the eigenvalues 
w| are real and positive except for the first several steps 
of the iteration. By taking U^{f7£,\£,) and J7'(cr^|<^) with 
^ € 1, • ■ • , M, we construct the projection operators 
and UK 

The renormalization procedure is described as follows. 
On the basis of the variational principle for the partition 
function,*^ we use the operator to renormalize the 
spin variables on each "edge" of X^ , and the operator 
is used for X'^. Thus the renormalized tensors X^ 
and X^ are obtained as 



X' 



X' 



^ij ^i'j ^i'j' ^ij' 
^ij ^i'j ^i'j' ^ij' 



(23) 




fig. 6 N. Maeshima 



Fig. 6. The vertical density matrix of the 3D classical model 
mapped from the 2D XXZ model. Open circles (rectangles) rep- 
resent the fixed spin variables in eq. (21), and filled circles (rect- 
angles) represent the summed variables. 



and 



X' 



E 



w 



X X^ {Cij^i'j£.i'j'^2j')U''{ai.iCij\^ij)U''{cri<j^i'j\^i'j) 



X ^/'■(CTi'j'C^'i'IC*'j')C^''(c^y'C»j'l6i')- 



(24) 



The transposed tensor is obtained by using the relation 
eq. (20). 

Here, it should be noted that the forms of the tensors 
alternate by one renormalization procedure; the new ten- 
sor X^ has only the auxiliary variables as the indices, 
whereas X^ has both the bare spins and auxiliary vari- 
ables. Therefore, X^^^^ has the same form as X-^^'^\ 
which suggests that we should interchange the tensors as 



new 1 



and 



new 1 



(25) 



to obtain new tensors X^^^ and X^^^ for the next step. 

The method of calculating the VDM p is essentially 
the same as that for the 3D Ising model. Define two 
tensors 



-x^ {v.^)x^ 



and 



(26) 



(27) 
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and regard and as the Boltzmann weights of a 
new 2D classical statistical system. Then we can calcu- 
late the VDM and various observables using the stan- 
dard DMRG method for 2D classical statistical sys- 
tems''^ or the corner transfer matrix renormalization 
group (CTMRG).''^ It should be also noted that, in the 
DMRG for the 2D classical system (or in the CTMRG), 
the number of the retained basis is denoted by m, which 
is discriminated from the number of state M of the aux- 
iliary spin variables ^ and rj. 

Finally some technical details of the VDMA are men- 
tioned. Since the VDMA deals with infinite-size systems, 
we can calculate the spontaneous staggered magnetiza- 
tion directly by applying the finite staggered magnetic 
field. In actual calculation the staggered magnetic field 
is introduced as follows. Only at the initial step of the 
VDMA we use the Boltzmann weight defined as 



fig.7 N.Maeshima 



initial 



(cTyl exp[e/is(y)(i3's)]|<7ij), (28) 



where the local Hamiltonian is 



h 



Hg {SI 



qZ 



TTX ( QX nx QX I nx \ fOQ\ 



and = {H^,Hy = 0, iJf ) is the staggered field. 

How the staggered field is applied influences converged 
results and the speed of the convergence. We should 
not apply the strong staggered field that conflicts the 
anisotropy of the system. For example, if we set >> 1 
and = for the Ising-like model {g < 1), the conver- 
gence often becomes very slow because such an inappro- 
priate field can lead the result to a quasi-stable XY-Neel 
state, not to the Ising-like Neel ground state. 

The author also gives a comment on the initial ten- 
sor w{t). It is found that w{t) does not have a strong 
influence to the converged results and the speed of the 
convergence if the initial TPS made from w{t) is not or- 
thogonal to IV'max)- The only one point of notice is that 
we must avoid such u'(t) as eflfectively apply a strong 
staggered field that conflicts the anisotropy. 

3. Results 

3. 1 The unfrustrated region (jj. = 0) 

First, the VDMA calculation has been performed for 
the 2D XXZ model on the square lattice {jd = 0) . This 
model has been studied intensively with several numeri- 
cal methods and thus reliable data are available. Hence 
we here compare the VDMA results with those already- 
known data to investigate the numerical accuracy of the 
VDMA. 

Figure 7 shows the e-dependence of the ground state 
energy per plaquette E and the staggered magnetiza- 
tion per site Mg for (M, m)=(2, 8) (open circles) and 
(3,9) (crosses). Here is the z-component {M^) of the 
staggered magnetization for the Ising-like model g < I, 
the x-component (Mf ) for the XY-like model g > 1, or 
the absolute value {^/{M^y + (M|)2) for the Heisen- 
berg model g = 1. It has been confirmed that the con- 
vergence with respect to m is sufiicient in the whole g 
region. 





d 0.1 6.2 




b ' ' 0.1 ' ' 0.2 ■ 



■ ■ 0.1' ■ ■ 0.2' 




■ ■ 0.1' . ■ 0.2' 



Fig. 7. Trotter extrapolation of the energy E and tlie staggered 
magnetization Mg. Open circles denote results for (M, m)=(2, 8), 
and crosses for (3,9). The solid lines are fits of data. 



Table I. Results of the ground state energy E and the staggered 



maj 


rnetization Ms- 


E 




g 


(M, m)=(2,8) 


(M,m)-(3,9) 


QMC 


0.4 


-1.0532 


-1.0532 


-1.056 '-'^ 


0.8 


-1.2106 


-1.2110 


-1.214 25 


1.0 


-1.3250 


-1.3276 


-1.33887 27 


1.2 


-1.5314 


-1.5329 


-1.538 25 


2.0 


-2.3871 


-2.3884 


-2.4 23 


Mg 


g 


(M,m)=(2,8) 


(M,m)=(3,9) 


QMC 


0.4 


0.482 


0.482 


0.48 25 


0.8 


0.432 


0.431 


0.41 25 


1.0 


0.396 


0.393 


0.3070 27 


2.0 


0.433 


0.433 


0.41 24 



When using numerical methods based on the Suzuki- 
Trotter decomposition, the extrapolation for e — > 
should be performed to obtain the true values of physical 
quantities. Then the following function 



0(e) = 0(e = 0) -h a2e2 + 046^ 



(30) 



where a2 and 04 are constants, is commonly employed 
for an observable O.^"^ The E-e^ plot always shows the 
nearly linear function (see Fig. 7 (a,c,e) ). Thus eq. (30) 
can be used for extrapolation of the energy in the whole g 
region. For the staggered magnetization Mg, in the Ising- 
like region {g < 0.6) Mg-e"^ plot also shows almost linear 
dependence. The g = 0.4 case is illustrated in Fig. 7 (b). 
By contrast, for the (nearly) critical systems {g > 0.6), 
the almost linear dependence cannot be observed, which 
is attributed to the insuflicient M-convergence for each 
e (see Fig. 7 (d),(f)). Then Ms at e = 0.1 are used as a 
rough estimate of Ms(0). 
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Fig. 8. The staggered magnetization for finite j^. 

In Table I, results of E and M, are compared with sev- 
eral QMC results. ^■'"^''''^^ The relative deviations from 
the QMC results are the order of 1%. For Mg, the quan- 
titative agreement with the QMC results is not so good 
as that of E; at the isotropic point, the VDMA re- 
sult of Ms is about 30% larger than a recent QMC re- 
sult. However, it has been confirmed that the VDMA 
results can reproduce the qualitative behavior of Mg, 
such as the phase transition between the Ising-like Neel 
phase {M^ ^ 0, = 0) and the XY-like Neel phase 
(M| = 0,M^ 7^ 0) at g = 1. 

Here, the author gives a comment on the calculation 
time of the VDMA. The dominant part of the VDMA 
calculation is the implementation of the DMRG (or 
CTMRG) for the 2D classical system, and the calcula- 
tion time of this part scales as x M*. Hence it is 
highly sensitive to M. For example, one VDMA calcula- 
tion with (M, m) = (3,9) takes dozens of hours of work- 
station time on an Alpha 21264 workstation, although 
that with (M, m) = (2, 8) takes only one or two hours. 
Thus it is unreasonable to further increase M. 

3.2 The frustrated region (jd =/= 0) 

In this subsection the VDMA results for the frustrated 
region {ja ^ 0) are presented. The most important topic 
in this region is the presence of a magnetically disordered 
phase at the "2D pyrochlore" point {jd = g = 
In addition, in the isotropic case {g = 1) the disordered 
phase is considered to survive up to jd = 0.75.'^"' ''^ In the 
following the chief object of interest is the system in the 
anisotropic region {g ^ 1), which has not been explored. 

In Fig. 8, the staggered magnetization Ms is plotted 
as a function of jd for the Ising-like system {g = 0.4), 
the isotropic system {g = 1), and the XY-like system 
{g = 1.4). VDMA results shown in this subsection are 
obtained by calculations with (M, m)=(2, 16). The ex- 
trapolation procedure for e — > is same as that of the 
unfrustrated case. The obtained results show that the fi- 
nite staggered magnetization Ms remains in the range 
< id < 1 for all g. For g = 1, this is contradict with 
the observation of the magnetically disordered phase for 
jd > 0.75.30.37 



Fig. 9. The ground state energy E. Symbols denote the VDMA 
results and the solid lines are guides to the eye. The dashed, 
dotted, and short-dashed lines show Eygc- 

The disordered phase observed at the 2D pyrochlore 
point is called the valence bond crystal (VBC),3^"35 
where the array of 4 spin singlets located on the void 
squares is the basic structure of the ground state. Here 
a simple argument is given based on the VBC picture 
to compare with the VDMA results. For this purpose, 
the VBC state is generalized for the anisotropic system. 
At the isotropic point the pure VBC state is a direct 
product of the ground state of a Heisenberg chain with 
only 4 spins. Hence for ^ 1 the generalized VBC state 
\'ipVBc) can be given as a direct product state of the 
ground states of the 4 spin chain with the XXZ-type 
interaction. Then the expectation value of the energy 
{EvBC = -1/4 - (v/l/4+ 2.g2)/2) can be obtained, 
which is an upper bound of the ground state energy. 
EvBC is independent of jd because of vanishing correla- 
tion of the next-nearest neighbor pairs interacting with 
jd- 

Figure 9 shows the VDMA results of the ground state 
energy E and the upper bound Evbc- It is foimd that 
the VDMA results for E is lower than EyBC in the re- 
gion where the frustration effect is not strong. However, 
the VDMA result for E exceeds Evbc around jd = 1, 
which implies that the VDMA calculation does not reach 
the true ground state there. This is the reason why the 
VDMA gives finite Ms even near jd = 1 for g = 1. 

The comparison between the VDMA results and the 
upper bound Ey bc leads to the phase diagram depicted 
in Fig. 10. Only from the VDMA, it is concluded that the 
whole parameter region is in the Neel phase. The g = 1 
line is the phase boundary between Ising-like phase and 
the XY-like one. Although the presence of the Neel phase 
is valid in the weakly frustrated region, it is quite doubt- 
ful in the shaded region in Fig. 10, where the generalized 
VBC state gives lower energy than the VDMA result. 
The remaining problem is what is the realized phase in 
the shaded region. Naturally, the VBC phase is the prime 
candidate. In particular, it is likely that the VBC phase 
extends around the g = 1 line within 0.75 < jd < 1.30-37 
However, it is unclear whether the VBC phase is realized 
far away from the 2D pyrochlore point. Hence, at present. 
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Fig. 10. Phase diagram of the S= 1/2 XXZ model on the checker- 
board lattice. In the shaded region, the generalized VBC state 
gives lower energy than the VDMA result. 
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Fig. 11. The valence bond crystal on the checkerboard lattice. 
The bold squares show the four-spin singlets. 

we cannot reach a definite conclusion on the shaded re- 
gion. Thus, further theoretical investigation is needed to 
complete the ground state phase diagram of this model. 

Finally, let us consider how we can obtain a good ap- 
proximation of the true ground state in the shaded region 
of Fig. 10. One possible solution in the TPS formulation 
is to construct a new variational TPS that can represent 
the pure VBC state. The simplest one of such a state can 
be made as follows. As shown in rcf.,"^^ the pure VBC 
state has the staggered pattern of the 4 spin singlets on 
the void squares (see Fig. 11). Hence the void squares 
can be divided into two groups C and D; the group C 
corresponds the void squares that have the singlet states 
and group D corresponds those in which 4 spins are un- 
corrclatcd. This structure of the VBC needs two types 
of tensors in the TPS representation; one constructs the 
4 spin singlet and another represents the uncorrelated 4 
spin state. Thus the TPS form of the pure VBC state is 
denoted as 

{\a\\VBC)= n ^cicTii) n ^^(^*^) (31) 

s{ij)eC s{ij)eD 



and the two types of tensors are given as 



-^c(o'l, 0'2, (73, 0-4) = < 



1 for (cti, (T2, 0-3, 0-4) 
= (Titi) or (itit) 

-i for (0-1,(72,(73,(74) 

= (Tni),(nir) (32) 
(iiTT),or am) 



for other 

(ai, (72, (73, 0-4), 

and 

-^_d((7i,CT2,o-3,(74) = 1 for any (171,(72,(73,(74). (33) 

Then a straightforward variational scheme based on the 
above argument is to take the several (or all) components 
of Xc and Xjj as variational parameters and to find the 
minimum of the variational energy with varying those 
parameters. 

It should be noted that such a variational TPS is the 
simplest one to represent the VBC state. The tensors Xc 
and Xd have only the bare spin variables cr, which cor- 
responds to the case of M = 1. Hence more complicated 
TPS having auxiliary variables with M > 2 would give 
lower Miriational energy because it has more variational 
parameters. 

The above discussion also tells us the reason why we 
have failed to reach the VBC state. In this study, all 
the tensors are put on the plaquettes that has the cross 
bond, not on the void plaquette. Hence the TPS used 
here docs not have the natural form to represent the 
pure VBC state. Thus, to complete the phase diagram, 
we probably need to use the TPS that reflects the spatial 
distribution of the 4 spin singlets in the VBC state. 

4. Conclusion 

In this paper the vertical density matrix algorithm 
(VDMA) is applied to the spin 1/2 XXZ model on the 
checkerboard lattice. The VDMA for this model needs 
several modifications in its implementation procedure be- 
cause of the lattice structure of the 3D classical spin sys- 
tem mapped by the Suzuki- Trotter transformation. 

The ground state energy and the staggered magnetiza- 
tion are calculated for the unfrustrated model to investi- 
gate the numerical accuracy of the VDMA. The obtained 
data show that the ground state energy is quantitatively 
reliable; the maximum deviation from the QMC results 
is about 1% for the unfrustrated Heisenberg model. For 
the staggered magnetization the VDMA only reproduces 
the qualitative behavior of the 2D XXZ model. 

In the frustrated region, the VDMA results are com- 
pared with the simple calculation based on the VBC 
state, which is considered to be realized at the 2D py- 
rochlore point. It is concluded that, in the weakly frus- 
trated region, the system is in the Neel ordered phase. 
By contrast, in the strongly frustrated region where the 
VDMA does not reach the true ground state, the realized 
phase cannot be identified. Hence to complete the phase 
diagram of this model, we need more reliable data for the 
strongly frustrated region in the thermodynamic limit. 



8 J. Phys. Soc. Jpn. 



Full Paper 



Nobuya Maeshima 



For this purpose, further improvement of the VDMA or 
a novel tensor product formulation would be useful. 
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